Assignment 7

Setup

knitr::opts_chunk$set(echo = TRUE, comment = "#>", dpi = 300)

for (f in list.files(here::here("src"), pattern = "R$", full.names = TRUE)) {
  source(f)
}

library(rstan)
library(tidybayes)
library(magrittr)
library(tidyverse)

theme_set(theme_classic() + theme(strip.background = element_blank()))

options(mc.cores = 2)
rstan_options(auto_write = TRUE)

drowning <- aaltobda::drowning
factory <- aaltobda::factory

Assignment 7

1. Linear model: drowning data with Stan

The provided data drowning in the ‘aaltobda’ package contains the number of people who died from drowning each year in Finland 1980–2019. A statistician is going to fit a linear model with Gaussian residual model to these data using time as the predictor and number of drownings as the target variable. She has two objective questions:

What is the trend of the number of people drowning per year? (We would plot the histogram of the slope of the linear model.)
What is the prediction for the year 2020? (We would plot the histogram of the posterior predictive distribution for the number of people drowning at \(\tilde{x} = 2020\).)

Corresponding Stan code is provided in Listing 1. However, it is not entirely correct for the problem. First, there are three mistakes. Second, there are no priors defined for the parameters. In Stan, this corresponds to using uniform priors.

a) Find the three mistakes in the code and fix them. Report the original mistakes and your fixes clearly in your report. Include the full corrected Stan code in your report.

Declaration of sigma on line 10 should be real<lower=0>.
Missing semicolon at the end of line 16.
On line 19, the prediction on new data does not use the new data in xpred. This has been changed to real ypred = normal_rng(alpha + beta*xpred, sigma);.

Below is a copy of the final model. The full Stan file is at models/assignment07-drownings.stan.

data {
  int<lower=0> N;  // number of data points
  vector[N] x;     // observation year
  vector[N] y;     // observation number of drowned
  real xpred;      // prediction year
}
parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;  // fix: 'upper' should be 'lower'
}
transformed parameters {
  vector[N] mu = alpha + beta*x;
}
model {
  y ~ normal(mu, sigma);  // fix: missing semicolor
}
generated quantities {
  real ypred = normal_rng(alpha + beta*xpred, sigma);  // fix: use `xpred`
}

b) Determine a suitable weakly-informative prior \(\text{Normal}(0,\sigma_\beta)\) for the slope \(\beta\). It is very unlikely that the mean number of drownings changes more than 50 % in one year. The approximate historical mean yearly number of drownings is 138. Hence, set \(\sigma_\beta\) so that the following holds for the prior probability for \(\beta\): \(Pr(−69 < \beta < 69) = 0.99\). Determine suitable value for \(\sigma_\beta\) and report the approximate numerical value for it.

x <- rnorm(1e5, 0, 26)
print(mean(-69 < x & x < 69))

#> [1] 0.99205

plot_single_hist(x, alpha = 0.5, color = "black") + geom_vline(xintercept = c(-69, 69)) + labs(x = "beta")

c) Using the obtained σβ, add the desired prior in the Stan code.

From some trial and error, it seems that a prior of \(\text{Normal}(0, 26)\) should work. I have added this prior distribution to beta in the model at line 17.

beta ~ normal(0, 26);   // prior on `beta`

d) In a similar way, add a weakly informative prior for the intercept alpha and explain how you chose the prior.

To use the year directly as the values for \(x\) would lead to a massive value of \(\alpha\) because the values for \(x\) range from 1980 to 2019. Thus, it would be advisable to first center the year, meaning at the prior distribution for \(\alpha\) can be centered around the average of the number of drownings per year and a standard deviation near that of the actual number of drownings.

head(drowning)

#>   year drownings
#> 1 1980       149
#> 2 1981       127
#> 3 1982       139
#> 4 1983       141
#> 5 1984       122
#> 6 1985       120

print(mean(drowning$drownings))

#> [1] 134.35

print(sd(drowning$drownings))

#> [1] 28.48441

Therefore, I add the prior \(\text{Normal}(135, 50)\) to \(\alpha\) on line 16.

alpha ~ normal(135, 50);   // prior on `alpha`

data <- list(
  N = nrow(drowning),
  x = drowning$year - mean(drowning$year),
  y = drowning$drownings,
  xpred = 2020 - mean(drowning$year)
)
drowning_model <- stan(
  here::here("models", "assignment07-drownings.stan"),
  data = data
)

variable_post <- spread_draws(drowning_model, alpha, beta) %>%
  pivot_longer(c(alpha, beta), names_to = "variable", values_to = "value")
head(variable_post)

#> # A tibble: 6 × 5
#>   .chain .iteration .draw variable  value
#>    <int>      <int> <int> <chr>     <dbl>
#> 1      1          1     1 alpha    131.
#> 2      1          1     1 beta      -1.30
#> 3      1          2     2 alpha    138.
#> 4      1          2     2 beta      -1.27
#> 5      1          3     3 alpha    131.
#> 6      1          3     3 beta      -1.05

variable_post %>%
  ggplot(aes(x = .iteration, y = value, color = factor(.chain))) +
  facet_grid(rows = vars(variable), scales = "free_y") +
  geom_path(alpha = 0.5) +
  scale_x_continuous(expand = expansion(c(0, 0))) +
  scale_y_continuous(expand = expansion(c(0.02, 0.02))) +
  labs(x = "iteration", y = "value", color = "chain")

variable_post %>%
  ggplot(aes(x = value)) +
  facet_grid(cols = vars(variable), scales = "free_x") +
  geom_histogram(color = "black", alpha = 0.3, bins = 30) +
  scale_x_continuous(expand = expansion(c(0.02, 0.02))) +
  scale_y_continuous(expand = expansion(c(0, 0.02)))

spread_draws(drowning_model, ypred) %$%
  plot_single_hist(ypred, alpha = 0.3, color = "black") +
  labs(x = "predicted number of drownings in 2020")

red <- "#C34E51"

bayestestR::describe_posterior(drowning_model, ci = 0.89, test = c()) %>%
  as_tibble() %>%
  filter(str_detect(Parameter, "mu")) %>%
  select(Parameter, Median, CI_low, CI_high) %>%
  janitor::clean_names() %>%
  mutate(idx = row_number()) %>%
  left_join(drowning %>% mutate(idx = row_number()), by = "idx") %>%
  ggplot(aes(x = year)) +
  geom_point(aes(y = drownings), data = drowning, color = "#4C71B0") +
  geom_line(aes(y = median), color = red, size = 1.2) +
  geom_smooth(
    aes(y = ci_low),
    method = "loess",
    formula = "y ~ x",
    linetype = 2,
    se = FALSE,
    color = red,
    size = 1
  ) +
  geom_smooth(
    aes(y = ci_high),
    method = "loess",
    formula = "y ~ x",
    linetype = 2,
    se = FALSE,
    color = red,
    size = 1
  ) +
  labs(x = "year", y = "number of drownings (mean ± 89% CI)")

2. Hierarchical model: factory data with Stan

The factory data in the ‘aaltobda’ package contains quality control measurements from 6 machines in a factory (units of the measurements are irrelevant here). In the data file, each column contains the measurements for a single machine. Quality control measurements are expensive and time-consuming, so only 5 measurements were done for each machine. In addition to the existing machines, we are interested in the quality of another machine (the seventh machine).

For this problem, you’ll use the following Gaussian models:

a separate model, in which each machine has its own model
a pooled model, in which all measurements are combined and there is no distinction between machines
a hierarchical model, which has a hierarchical structure as described in BDA3 Section 11.6

As in the model described in the book, use the same measurement standard deviation \(\sigma\) for all the groups in the hierarchical model. In the separate model, however, use separate measurement standard deviation \(\sigma_j\) for each group \(j\). You should use weakly informative priors for all your models.

Complete the following questions for each of the three models (separate, pooled, hierarchical).

a) Describe the model with mathematical notation. Also describe in words the difference between the three models.

Separate model: The separate model is described below where each machine has its own centrality \(\mu\) and dispersion \(\sigma\) parameters that do not influence the parameters of the other machines.

\[ y_{ij} \sim N(\mu_j, \sigma_j) \\ \mu_j \sim N(0, 1) \\ \sigma_j \sim \text{Inv-}\chi^2(10) \]

Pooled model: The pool model is described below where there is no distinction between the models but instead a single set of parameters for all of the data.

\[ y_{i} \sim N(\mu, \sigma) \\ \mu \sim N(0, 1) \\ \sigma \sim \text{Inv-}\chi^2(10) \]

Hierarchical model: The hierarchical model is described below where each machine has its own centrality \(\mu\) parameter which are linked through a hyper-prior distribution from which they are drawn. The machines will all share a common dispersion paramete \(\sigma\)

\[ y_{ij} \sim N(\mu_j, \sigma_j) \\ \mu_j \sim N(\alpha, \tau) \\ \alpha \sim N(0, 1) \\ \tau \sim \text{HalfNormal}(2.5) \\ \sigma \sim \text{Inv-}\chi^2(10) \]

The separate model is effectively building a different linear model for each machine where as the pooled model treats all the measurements as coming from the same model. The hierarchical model is treating the machines as having come from a single, shared distribution.

b) Implement the model in Stan and include the code in the report. Use weakly informative priors for all your models.

print_model_code <- function(path) {
  for (l in readLines(path)) {
    cat(l, "\n")
  }
}

Separate model

separate_model_code <- here::here(
  "models", "assignment07_factories_separate.stan"
)
print_model_code(separate_model_code)

#> data {
#>   int<lower=0> N;  // number of data points per machine
#>   int<lower=0> J;  // number of machines
#>   vector[J] y[N];  // quality control data points
#> }
#>
#> parameters {
#>   vector[J] mu;
#>   vector<lower=0>[J] sigma;
#> }
#>
#> model {
#>   // priors
#>   for (j in 1:J) {
#>     mu[j] ~ normal(100, 10);
#>     sigma[j] ~ inv_chi_square(5);
#>   }
#>
#>   // likelihood
#>   for (j in 1:J){
#>     y[,j] ~ normal(mu[j], sigma[j]);
#>   }
#> }
#>
#> generated quantities {
#>   // Compute the predictive distribution for the sixth machine.
#>   real y6pred;
#>   y6pred = normal_rng(mu[6], sigma[6]);
#> }

separate_model_data <- list(
  y = factory,
  N = nrow(factory),
  J = ncol(factory)
)
separate_model <- rstan::stan(
  separate_model_code,
  data = separate_model_data,
  verbose = FALSE,
  refresh = 0
)
knitr::kable(
  bayestestR::describe_posterior(separate_model, ci = 0.89, test = NULL),
  digits = 3
)

Parameter	Median	CI	CI_low	CI_high	ESS	Rhat
mu[1]	85.044	0.89	74.117	96.354	3178.483	1.000
mu[2]	105.195	0.89	98.268	111.956	3612.836	1.000
mu[3]	90.033	0.89	83.033	97.666	3763.803	1.000
mu[4]	110.683	0.89	106.223	115.422	3386.517	1.000
mu[5]	91.381	0.89	85.301	97.883	3791.011	1.001
mu[6]	90.978	0.89	82.003	101.176	4415.799	1.000
y6pred	90.849	0.89	61.153	119.877	4083.663	1.000

Pooled model

pooled_model_code <- here::here("models", "assignment07_factories_pooled.stan")
print_model_code(pooled_model_code)

#> data {
#>   int<lower=0> N;  // number of data points
#>   vector[N] y;     // machine quality control data
#> }
#>
#> parameters {
#>   real mu;
#>   real<lower=0> sigma;
#> }
#>
#> model {
#>   // priors
#>   mu ~ normal(100, 10);
#>   sigma ~ inv_chi_square(5);
#>
#>   // likelihood
#>   y ~ normal(mu, sigma);
#> }
#>
#> generated quantities {
#>   real ypred;
#>   ypred = normal_rng(mu, sigma);
#> }

pooled_model_data <- list(
  y = unname(unlist(factory)),
  N = length(unlist(factory))
)
pooled_model <- rstan::stan(
  pooled_model_code,
  data = pooled_model_data,
  verbose = FALSE,
  refresh = 0
)

knitr::kable(
  bayestestR::describe_posterior(pooled_model, ci = 0.89, test = NULL),
  digits = 3
)

Parameter	Median	CI	CI_low	CI_high	ESS	Rhat
mu	93.585	0.89	88.888	99.123	2746.458	1.002
ypred	93.365	0.89	64.689	122.717	4146.168	1.000

Hierarchical model

hierarchical_model_code <- here::here(
  "models", "assignment07_factories_hierarchical.stan"
)
print_model_code(hierarchical_model_code)

#> data {
#>   int<lower=0> N;  // number of data points per machine
#>   int<lower=0> J;  // number of machines
#>   vector[J] y[N];  // quality control data points
#> }
#>
#> parameters {
#>   vector[J] mu;
#>   real<lower=0> sigma;
#>   real alpha;
#>   real<lower=0> tau;
#> }
#>
#> model {
#>   // hyper-priors
#>   alpha ~ normal(100, 10);
#>   tau ~ normal(0, 10);
#>
#>   // priors
#>   mu ~ normal(alpha, tau);
#>   sigma ~ inv_chi_square(5);
#>
#>   // likelihood
#>   for (j in 1:J){
#>     y[,j] ~ normal(mu[j], sigma);
#>   }
#> }
#>
#> generated quantities {
#>   // Compute the predictive distribution for the sixth machine.
#>   real y6pred;
#>   real mu7pred;
#>   y6pred = normal_rng(mu[6], sigma);
#>   mu7pred = normal_rng(alpha, tau);
#> }

hierarchical_model_data <- list(
  y = factory,
  N = nrow(factory),
  J = ncol(factory)
)
hierarchical_model <- rstan::stan(
  hierarchical_model_code,
  data = hierarchical_model_data,
  verbose = FALSE,
  refresh = 0
)
knitr::kable(
  bayestestR::describe_posterior(hierarchical_model, ci = 0.89, test = NULL),
  digits = 3
)

	Parameter	Median	CI	CI_low	CI_high	ESS	Rhat
2	mu[1]	81.424	0.89	71.244	91.954	1267.017	1.004
3	mu[2]	102.460	0.89	92.868	111.689	2062.469	1.001
4	mu[3]	89.675	0.89	81.240	99.006	3329.720	1.000
5	mu[4]	106.319	0.89	96.998	117.021	1257.313	1.003
6	mu[5]	91.321	0.89	82.252	100.030	3422.894	1.001
7	mu[6]	88.679	0.89	79.229	96.759	2560.124	1.002
1	alpha	94.370	0.89	87.169	101.616	2923.110	1.001
9	tau	10.501	0.89	4.634	17.616	984.551	1.004
10	y6pred	88.701	0.89	63.541	112.643	3800.454	1.001
8	mu7pred	94.159	0.89	75.719	114.333	3910.469	1.001

c) Using the model (with weakly informative priors) report, comment on and, if applicable, plot histograms for the following distributions:

the posterior distribution of the mean of the quality measurements of the sixth machine.
the predictive distribution for another quality measurement of the sixth machine.
the posterior distribution of the mean of the quality measurements of the seventh machine.

plot_hist_mean_of_sixth <- function(vals) {
  plot_single_hist(vals, bins = 30, color = "black", alpha = 0.3) +
    labs(x = "mean of 6th machine", y = "posterior density")
}

plot_hist_sixth_predictions <- function(vals) {
  plot_single_hist(vals, bins = 30, color = "black", alpha = 0.3) +
    labs(x = "posterior predictions for 6th machine", y = "posterior density")
}

plot_hist_mean_of_seventh <- function(vals) {
  plot_single_hist(vals, bins = 30, color = "black", alpha = 0.3) +
    labs(x = "mean of 7thth machine", y = "posterior density")
}

Separate model

plot_hist_mean_of_sixth(rstan::extract(separate_model)$mu[, 6])

plot_hist_sixth_predictions(rstan::extract(separate_model)$y6pred)

It is not possible to estimate the posterior for the mean of some new 7th machine because all machines are treated separately.

Pooled model

plot_hist_mean_of_sixth(rstan::extract(pooled_model)$mu)

plot_hist_sixth_predictions(rstan::extract(pooled_model)$ypred)

The predicted mean for a new machine is the same as the pooled mean \(mu\).

plot_hist_mean_of_seventh(rstan::extract(pooled_model)$mu)

Hierarchical model

plot_hist_mean_of_sixth(rstan::extract(hierarchical_model)$mu[, 6])

plot_hist_sixth_predictions(rstan::extract(hierarchical_model)$y6pred)

plot_hist_mean_of_seventh(rstan::extract(hierarchical_model)$mu7pred)

d) Report the posterior expectation for \(\mu_1\) with a 90% credible interval but using a \(\text{Normal}(0,10)\) prior for the \(\mu\) parameter(s) and a \(\text{Gamma}(1,1)\) prior for the \(\sigma\) parameter(s). For the hierarchical model, use the \(\text{Normal}(0, 10)\) and \(\text{Gamma}(1, 1)\) as hyper-priors.

(I’m going to skip this one, but come back to it if it is needed for future assignments.)

sessionInfo()

#> R version 4.1.1 (2021-08-10)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Big Sur 10.16
#>
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> attached base packages:
#> [1] stats     graphics  grDevices datasets  utils     methods
#> [7] base
#>
#> other attached packages:
#>  [1] forcats_0.5.1        stringr_1.4.0        dplyr_1.0.7
#>  [4] purrr_0.3.4          readr_2.0.1          tidyr_1.1.3
#>  [7] tibble_3.1.3         tidyverse_1.3.1      magrittr_2.0.1
#> [10] tidybayes_3.0.1      rstan_2.21.2         ggplot2_3.3.5
#> [13] StanHeaders_2.21.0-7
#>
#> loaded via a namespace (and not attached):
#>  [1] nlme_3.1-152         fs_1.5.0             matrixStats_0.61.0
#>  [4] lubridate_1.7.10     insight_0.14.4       httr_1.4.2
#>  [7] rprojroot_2.0.2      tensorA_0.36.2       tools_4.1.1
#> [10] backports_1.2.1      bslib_0.2.5.1        utf8_1.2.2
#> [13] R6_2.5.0             mgcv_1.8-36          DBI_1.1.1
#> [16] colorspace_2.0-2     ggdist_3.0.0         withr_2.4.2
#> [19] tidyselect_1.1.1     gridExtra_2.3        prettyunits_1.1.1
#> [22] processx_3.5.2       downlit_0.2.1        curl_4.3.2
#> [25] compiler_4.1.1       rvest_1.0.1          cli_3.0.1
#> [28] arrayhelpers_1.1-0   xml2_1.3.2           bayestestR_0.11.0
#> [31] labeling_0.4.2       posterior_1.1.0      sass_0.4.0
#> [34] scales_1.1.1         checkmate_2.0.0      aaltobda_0.3.1
#> [37] callr_3.7.0          digest_0.6.27        rmarkdown_2.10
#> [40] pkgconfig_2.0.3      htmltools_0.5.1.1    highr_0.9
#> [43] dbplyr_2.1.1         rlang_0.4.11         readxl_1.3.1
#> [46] rstudioapi_0.13      jquerylib_0.1.4      farver_2.1.0
#> [49] generics_0.1.0       svUnit_1.0.6         jsonlite_1.7.2
#> [52] distill_1.2          distributional_0.2.2 inline_0.3.19
#> [55] loo_2.4.1            Matrix_1.3-4         Rcpp_1.0.7
#> [58] munsell_0.5.0        fansi_0.5.0          abind_1.4-5
#> [61] lifecycle_1.0.0      stringi_1.7.3        yaml_2.2.1
#> [64] snakecase_0.11.0     pkgbuild_1.2.0       grid_4.1.1
#> [67] parallel_4.1.1       crayon_1.4.1         lattice_0.20-44
#> [70] splines_4.1.1        haven_2.4.3          hms_1.1.0
#> [73] knitr_1.33           ps_1.6.0             pillar_1.6.2
#> [76] codetools_0.2-18     clisymbols_1.2.0     stats4_4.1.1
#> [79] reprex_2.0.1         glue_1.4.2           evaluate_0.14
#> [82] V8_3.4.2             renv_0.14.0          RcppParallel_5.1.4
#> [85] modelr_0.1.8         vctrs_0.3.8          tzdb_0.1.2
#> [88] cellranger_1.1.0     gtable_0.3.0         datawizard_0.2.1
#> [91] assertthat_0.2.1     xfun_0.25            janitor_2.1.0
#> [94] broom_0.7.9          coda_0.19-4          ellipsis_0.3.2
#> [97] here_1.0.1